For this analysis I use the HCC data set with the relative water levels (counterclimatic and observed). Data can be accessed here

Global and Regional Water level changes (fig.13)

In this part we will make a plot of the regional observed water level changes and compare it to the global average. What we use for this is the hourly hcc waterlevel data that can be downloaded for every year seperately. In order to get a first idea of this data, we can take one example year: For this we extract the relevant variables and then take Stockholm as location.

setwd("/Users/hannascheiwe/Desktop/CLEWS2.0/uni_kram/CC_Adaptation/data/hourly_data_relative_waterlevels_observed")

wl1997 <- nc_open("hcc_obsclim_waterlevel_global_hourly_1997.nc")
print(wl1997)
## File hcc_obsclim_waterlevel_global_hourly_1997.nc (NC_FORMAT_NETCDF4):
## 
##      3 variables (excluding dimension variables):
##         double waterlevel[stations,time]   (Contiguous storage)  
##             comment: This variable represents the relative coastal water level as a combination of reconstructed geocentric water level data based on Dangendorf et al. (2019), vertical land motion (VLM) based on Oelsmann et al. (2023) and hourly storm-tide variations based on Muis et al. (2020). The data is aligned with the geocentric water level such that the last time step in the record is identical. It therefore describes the water level relative to a virtual tide gauge that is at the height of the geocentric water level at the last date in the record.
##             units: mm
##             positive: up
##             coordinates: lat lon
##             long_name: Relative Water Level
##         double lat[stations]   (Contiguous storage)  
##             standard_name: latitude
##             long_name: Latitude
##             units: degrees_north
##             axis: Y
##         double lon[stations]   (Contiguous storage)  
##             standard_name: longitude
##             long_name: Longitude
##             units: degrees_east
##             axis: X
## 
##      2 dimensions:
##         stations  Size:18719 (no dimvar)
##         time  Size:8760 
##             long_name: Time
##             axis: T
##             units: hours since 1979-01-01
##             calendar: proleptic_gregorian
## 
##     9 global attributes:
##         Convention: CF-1.9
##         featureType: timeSeries
##         references: The dataset is described in Treu et al. (2023) <https://doi.org/10.5194/essd-2023-112>. The source datasets are described in Dangendorf et al. (2019) <https://doi.org/10.1038/s41558-019-0531-8>, Muis et al. (2020) <https://doi.org/10.3389/fmars.2020.00263> and Oelsmann et al. (2023) <https://doi.org/10.5281/zenodo.8308347>
##         institution: Potsdam Institute for Climate Impact Research (PIK)
##         project: Inter-Sectoral Impact Model Intercomparison Project phase 3a (ISIMIP3a)
##         contact: ISIMIP cross-sectoral science team <info@isimip.org> <https://www.isimip.org>
##         isimip_id: d193c353-beee-4e8e-bdd4-84d7a6b74297
##         grid_info: The irregular coastal grid points are distributed equidistantly along a smoothed global coastlines with a distance of 50 km globally and 10km in Europe with additional points added for tide gauge stations (described in Muis et all. 2019, https://doi.org/10.3389/fmars.2020.00263).
##         title: Hourly coastal relative water levels
water_level_1997 <- ncvar_get(wl1997, "waterlevel")
lat <- ncvar_get(wl1997, "lat")
lon <- ncvar_get(wl1997, "lon")
time_values_1997 <- ncvar_get(wl1997, "time")


stockholm_lat <- 59.3293
stockholm_lon <- 18.0686

distances_1997 <- sqrt((lat - stockholm_lat)^2 + (lon - stockholm_lon)^2)
closest_station_idx_1997 <- which.min(distances_1997)

stockholm_water_level_1997 <- water_level_1997[closest_station_idx_1997 , ]  


# Calculate global average water level for each time step (average across all stations)
global_average_water_level_1997 <- rowMeans(water_level_1997, na.rm = TRUE)

mean_water_level_all_stations_1997 <- apply(water_level_1997, 2, mean, na.rm = TRUE) 
#--> calculates the average water level across all stations for each time step


global_avg_sea_level_1997 <- mean(mean_water_level_all_stations_1997, na.rm = TRUE)

Okay now we do have the average values we were looking for, but the goal is to have a time series over the whole period. So the next step will be to get the mean value of all stations for the respective year in order to get the development of the global average sea level change. For this we can use a loop that applies all these steps for the different data. Here it is important that all the files are stored in the same folder.

global_avg_sea_levels <- vector("numeric", length = 19) 

stockholm_water_levels <- vector("list", length = 19)  

for (year in 1997:2015) {
  file_name <- paste0("hcc_obsclim_waterlevel_global_hourly_", year, ".nc")
  if (file.exists(file_name)) {
    dat <- nc_open(file_name)
    water_level <- ncvar_get(dat, "waterlevel")
    lat <- ncvar_get(dat, "lat")
    lon <- ncvar_get(dat, "lon")
    time_values <- ncvar_get(dat, "time")
    distances <- sqrt((lat - stockholm_lat)^2 + (lon - stockholm_lon)^2)
    closest_station_idx <- which.min(distances)
    stockholm_water_level <- water_level[closest_station_idx , ]
    stockholm_water_levels[[year - 1996]] <- stockholm_water_level
    global_average_water_level <- rowMeans(water_level, na.rm = TRUE)
    mean_water_level_all_stations <- apply(water_level, 2, mean, na.rm = TRUE)
    global_avg_sea_level <- mean(mean_water_level_all_stations, na.rm = TRUE)
    

    global_avg_sea_levels[year - 1996] <- global_avg_sea_level 
  } else {
    
    cat("File for year", year, "not found. Skipping...\n")
  }
}
## File for year 2001 not found. Skipping...
print(global_avg_sea_levels)
##  [1] -9.452900  3.475170  9.975784  9.191684  0.000000  2.850987  7.044084
##  [8]  9.830954  8.804978 11.528670 19.242645 22.658903 19.638729 27.157361
## [15] 36.809982 33.308854 33.948932 33.369887 40.109891

Stockholm

There might be some repetition in the code and we could probably make it way easier but this is how it worked out in the end, but feel free to play around with it. Also we will iteratively plot the time series of the region to the global average.

years <- 1997:2015
yearly_stockholm_water_levels <- vector("numeric", length = 19)

for (year in 1997:2015) {
  file_name <- paste0("hcc_obsclim_waterlevel_global_hourly_", year, ".nc")
  if (file.exists(file_name)) {
    dat <- nc_open(file_name)
    water_level <- ncvar_get(dat, "waterlevel")
    lat <- ncvar_get(dat, "lat")
    lon <- ncvar_get(dat, "lon")

    distances <- sqrt((lat - stockholm_lat)^2 + (lon - stockholm_lon)^2)
    closest_station_idx <- which.min(distances)
    
    
    stockholm_water_level <- water_level[closest_station_idx , ]
    
   
    yearly_stockholm_water_levels[year - 1996] <- mean(stockholm_water_level, na.rm = TRUE)
    
    
    global_average_water_level <- rowMeans(water_level, na.rm = TRUE)
    
    
    mean_water_level_all_stations <- apply(water_level, 2, mean, na.rm = TRUE)
    
    
    global_avg_sea_levels[year - 1996] <- mean(mean_water_level_all_stations, na.rm = TRUE)  
  } else {
    cat("File for year", year, "not found. Skipping...\n")
  }
}
## File for year 2001 not found. Skipping...

Manila

manila_lat <- 14.5995  
manila_lon <- 120.9842 

yearly_manila_water_levels <- vector("numeric", length = 19)


for (year in 1997:2015) {
  file_name <- paste0("hcc_obsclim_waterlevel_global_hourly_", year, ".nc")
  if (file.exists(file_name)) {
    dat <- nc_open(file_name)
    water_level <- ncvar_get(dat, "waterlevel")
    lat <- ncvar_get(dat, "lat")
    lon <- ncvar_get(dat, "lon")
    distances <- sqrt((lat - manila_lat)^2 + (lon - manila_lon)^2)
    closest_station_idx <- which.min(distances)
    manila_water_level <- water_level[closest_station_idx , ]
    yearly_manila_water_levels[year - 1996] <- mean(manila_water_level, na.rm = TRUE)
    global_average_water_level <- rowMeans(water_level, na.rm = TRUE)
    
    mean_water_level_all_stations <- apply(water_level, 2, mean, na.rm = TRUE)
    

    global_avg_sea_levels[year - 1996] <- mean(mean_water_level_all_stations, na.rm = TRUE) 
  } else {

    cat("File for year", year, "not found. Skipping...\n")
  }
}
## File for year 2001 not found. Skipping...
######add manila to the plot

data <- data.frame(
  Year = rep(years, 3),
  Water_Level = c(global_avg_sea_levels, yearly_stockholm_water_levels, yearly_manila_water_levels),
  Series = rep(c("Global Average", "Stockholm", "Manila"), each = length(years))
)


ggplot(data, aes(x = Year, y = Water_Level, color = Series)) +
  geom_line(size = 1.2) +
  labs(
    title = "Global, Stockholm and Manila Water Level Changes (1997-2015)",
    x = "Year",
    y = "Water Level (mm)"
  ) +
  scale_color_manual(values = c("blue", "red", "orange")) +
  theme_minimal() +
  theme(legend.title = element_blank())

Charlottetown

charlottetown_lat <- 46.2382
charlottetown_lon <- -63.1311

yearly_charlottetown_water_levels <- vector("numeric", length = 19)


for (year in 1997:2015) {
  file_name <- paste0("hcc_obsclim_waterlevel_global_hourly_", year, ".nc")
  if (file.exists(file_name)) {
    dat <- nc_open(file_name)
    water_level <- ncvar_get(dat, "waterlevel")
    lat <- ncvar_get(dat, "lat")
    lon <- ncvar_get(dat, "lon")
    distances <- sqrt((lat - charlottetown_lat)^2 + (lon - charlottetown_lon)^2)
    closest_station_idx <- which.min(distances)
    charlottetown_water_level <- water_level[closest_station_idx , ]
    yearly_charlottetown_water_levels[year - 1996] <- mean(charlottetown_water_level, na.rm = TRUE)
    global_average_water_level <- rowMeans(water_level, na.rm = TRUE)
    mean_water_level_all_stations <- apply(water_level, 2, mean, na.rm = TRUE)
    

    global_avg_sea_levels[year - 1996] <- mean(mean_water_level_all_stations, na.rm = TRUE)  
  } else {
    cat("File for year", year, "not found. Skipping...\n")
  }
}
## File for year 2001 not found. Skipping...
# plot
data <- data.frame(
  Year = rep(years, 4),
  Water_Level = c(global_avg_sea_levels, yearly_stockholm_water_levels, yearly_manila_water_levels, yearly_charlottetown_water_levels),
  Series = rep(c("Global Average", "Stockholm", "Manila", "Charlottetown"), each = length(years))
)


ggplot(data, aes(x = Year, y = Water_Level, color = Series)) +
  geom_line(size = 1.2) +
  labs(
    title = "Global, Stockholm, Manila, and Charlottetown Water Level Changes (1997-2015)",
    x = "Year",
    y = "Water Level (mm)"
  ) +
  scale_color_manual(values = c("blue", "red", "orange", "green")) +
  theme_minimal() +
  theme(legend.title = element_blank())

final combination

data <- data.frame(
  year = rep(years, 4),
  water_level = c(global_avg_sea_levels, yearly_stockholm_water_levels, 
                  yearly_manila_water_levels, yearly_charlottetown_water_levels),
  region = rep(c("Global Average", "Stockholm", "Manila", "Charlottetown"), each = length(years))
)


data$region <- factor(data$region, levels = c("Global Average", "Stockholm", "Manila", "Charlottetown"))


water_level_changes_plot.png <-ggplot(data, aes(x = year, y = water_level, color = region, linetype = region)) +
  geom_line(size = 1.7) +  
  scale_color_manual(values = c("blue", "red", "orange", "green")) +  
  scale_linetype_manual(values = c("solid", "dashed", "dashed", "dashed")) +  
  labs(title = "Global and Regional Water Level Changes (1997-2015)", 
       subtitle = "Comparison of global average water level and regional locations",
       x = "Year", 
       y = "Water Level (mm)") +
  theme_minimal() + 
  theme(
    legend.position = "top",  
    legend.title = element_blank(),  
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10),
    plot.title = element_text(size = 16, face = "bold"),
    plot.subtitle = element_text(size = 12, face = "italic")
  ) +
  scale_x_continuous(breaks = seq(1997, 2015, by = 1)) +  
  scale_y_continuous(expand = c(0, 0), limits = c(min(data$water_level, na.rm = TRUE) - 10, 
                                                  max(data$water_level, na.rm = TRUE) + 10))  

print(water_level_changes_plot.png)

#ggsave("water_level_changes_plot.png", width = 10, height = 6, dpi = 300)

Observed and Counterclimatic Coastal Water Levels

Load the data and extract the water level information. Then calculate the mean of the waterlevels over time. This is done for the counterclimatic and observed values. By print() we can take a look at the data and gain some information on what is included, and most importantly we can see the names of the variables we would like to extract.

Load and Process Observed Water Levels (fig.14)

observed_hcc_relative <- nc_open("hcc_obsclim_waterlevel_global_monthly_1901_1978.nc")
print(observed_hcc_relative)
## File hcc_obsclim_waterlevel_global_monthly_1901_1978.nc (NC_FORMAT_NETCDF4):
## 
##      3 variables (excluding dimension variables):
##         double lat[stations]   (Contiguous storage)  
##             standard_name: latitude
##             long_name: Latitude
##             units: degrees_north
##             axis: Y
##         double lon[stations]   (Contiguous storage)  
##             standard_name: longitude
##             long_name: Longitude
##             units: degrees_east
##             axis: X
##         double waterlevel[stations,time]   (Contiguous storage)  
##             comment: This variable represents the monthly relative coastal water level as a combination of reconstructed geocentric water level data based on Dangendorf et al. (2019) and vertical land motion (VLM) based on Oelsmann et al. (2023). The data is aligned with the HCC hourly relative coastal water levels.
##             coordinates: lon lat
##             long_name: Relative Water Level
## 
##      2 dimensions:
##         stations  Size:18719 (no dimvar)
##         time  Size:936 
##             long_name: Time
##             axis: T
##             units: days since 1900-01-16
##             calendar: proleptic_gregorian
## 
##     10 global attributes:
##         Convention: CF-1.9
##         featureType: timeSeries
##         reference: https://doi.org/10.1038/s41558-019-0531-8
##         references: The dataset is described in Treu et al. (2023) <https://doi.org/10.5194/essd-2023-112>. The source datasets are described in Dangendorf et al. (2019) <https://doi.org/10.1038/s41558-019-0531-8>, Muis et al. (2020) <https://doi.org/10.3389/fmars.2020.00263> and Oelsmann et al. (2023) <https://doi.org/10.5281/zenodo.8308347>
##         institution: Potsdam Institute for Climate Impact Research (PIK)
##         project: Inter-Sectoral Impact Model Intercomparison Project phase 3a (ISIMIP3a)
##         contact: ISIMIP cross-sectoral science team <info@isimip.org> <https://www.isimip.org>
##         isimip_id: 79f4a736-4f0d-4b24-a6af-88d99310b200
##         grid_info: The irregular coastal grid points are distributed equidistantly along a smoothed global coastlines with a distance of 50 km globally and 10km in Europe with additional points added for tide gauge stations (described in Muis et all. 2019, https://doi.org/10.3389/fmars.2020.00263).
##         title: Monthly coastal relative water levels
rel_waterlevel_observed <- ncvar_get(observed_hcc_relative, "waterlevel")

mean_rel_obs <- rowMeans(rel_waterlevel_observed, na.rm = TRUE)

Load and Process Counterclimatic Water Levels

counterclim_hcc_rel <- nc_open("hcc_counterclim_waterlevel_global_monthly_1901_1978.nc")
print(counterclim_hcc_rel)
## File hcc_counterclim_waterlevel_global_monthly_1901_1978.nc (NC_FORMAT_NETCDF4):
## 
##      3 variables (excluding dimension variables):
##         double lat[stations]   (Contiguous storage)  
##             standard_name: latitude
##             long_name: Latitude
##             units: degrees_north
##             axis: Y
##         double lon[stations]   (Contiguous storage)  
##             standard_name: longitude
##             long_name: Longitude
##             units: degrees_east
##             axis: X
##         double waterlevel[stations,time]   (Contiguous storage)  
##             comment: This variable represents the counterfactual monthly relative coastal water level.It is derived by removing the quadratic trend since 1900 from the factual monthly relative coastal water level. The factual data is a combination of reconstructed geocentric water level data based on Dangendorf et al. (2019) and vertical land motion (VLM) based on Oelsmann et al. (2023).
##             coordinates: lon lat
##             long_name: Relative Water Level
## 
##      2 dimensions:
##         stations  Size:18719 (no dimvar)
##         time  Size:936 
##             long_name: Time
##             axis: T
##             units: days since 1900-01-16
##             calendar: proleptic_gregorian
## 
##     10 global attributes:
##         Convention: CF-1.9
##         featureType: timeSeries
##         reference: https://doi.org/10.1038/s41558-019-0531-8
##         references: The dataset is described in Treu et al. (2023) <https://doi.org/10.5194/essd-2023-112>. The source datasets are described in Dangendorf et al. (2019) <https://doi.org/10.1038/s41558-019-0531-8>, Muis et al. (2020) <https://doi.org/10.3389/fmars.2020.00263> and Oelsmann et al. (2023) <https://doi.org/10.5281/zenodo.8308347>
##         institution: Potsdam Institute for Climate Impact Research (PIK)
##         project: Inter-Sectoral Impact Model Intercomparison Project phase 3a (ISIMIP3a)
##         contact: ISIMIP cross-sectoral science team <info@isimip.org> <https://www.isimip.org>
##         isimip_id: 65151046-bbdd-436b-b4b1-2809e39eb9a1
##         grid_info: The irregular coastal grid points are distributed equidistantly along a smoothed global coastlines with a distance of 50 km globally and 10km in Europe with additional points added for tide gauge stations (described in Muis et all. 2019, https://doi.org/10.3389/fmars.2020.00263).
##         title: Counterfactual monthly coastal relative water levels
rel_waterlevel_counter <- ncvar_get(counterclim_hcc_rel, "waterlevel")

mean_rel_counter <- rowMeans(rel_waterlevel_counter, na.rm = TRUE)

Focus on Specific Locations

In the next steps, we look at three different locations (Stockholm, Charlottetown & Manila) and plot the different sea level trends that can be seen in these regions.

Map of Selected Locations

First, let´s create a map with the locations highlighted:

stockholm_lat <- 59.3293 
stockholm_lon <- 18.0686 

stockholm_location <- data.frame(
  lon = stockholm_lon,
  lat = stockholm_lat,
  name = "Stockholm"
)

charlottetown_lat <- 46.2382
charlottetown_lon <- -63.1311

charlottetown_location <- data.frame(
  lon = charlottetown_lon,
  lat = charlottetown_lat,
  name = "Charlottetown"
)

manila_lat <- 14.5995  
manila_lon <- 120.9842 

manila_location <- data.frame(
  lon = manila_lon,
  lat = manila_lat,
  name = "Manila"
)



world_plot <- ggplot() +
  geom_sf(data = world, fill = "gray90", color = "white") +
  geom_point(
    data = charlottetown_location, 
    aes(x = lon, y = lat), 
    color = "red", 
    size = 2
  ) +   
  geom_text(
    data = charlottetown_location, 
    aes(x = lon, y = lat, label = name), 
    hjust = -0.2, vjust = 0, size = 4, color = "red"
  ) +
  geom_point(
    data = stockholm_location, 
    aes(x = lon, y = lat), 
    color = "blue", 
    size = 2
  ) +
  geom_text(
    data = stockholm_location, 
    aes(x = lon, y = lat, label = name), 
    hjust = -0.2, vjust = 0, size = 4, color = "blue"
  ) + 
  geom_point(
    data = manila_location, 
    aes(x = lon, y = lat), 
    color = "purple", 
    size = 2
  ) +
  geom_text(
    data = manila_location, 
    aes(x = lon, y = lat, label = name), 
    hjust = -0.2, vjust = 0, size = 4, color = "purple"
  ) +
  labs(
    title = "Locations on the World Map",
    x = "Longitude",
    y = "Latitude"
  ) +
  coord_sf(xlim = c(-180, 180), ylim = c(-90, 90), expand = FALSE) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold"),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )

print(world_plot)

Next, we can plot time series for the single locations. This is done for the counterclimatic and observed data. On top of that we can calculate sea level trends for each data set and visualize this.